Symmetric hyperbolicity and consistent boundary conditions for second-order 

Einstein equations 
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We present two families of first-order in time and second-order in space formulations of the 
Einstein equations (variants of the Arnowitt-Deser-Misner formulation) that admit a complete set of 
characteristic variables and a conserved energy that can be expressed in terms of the characteristic 
variables. The associated constraint system is also symmetric hyperbolic in this sense, and all 
characteristic speeds are physical. We propose a family of constraint-preserving boundary conditions 
that is applicable if the boundary is smooth with tangential shift. We conjecture that the resulting 
initial-boundary value problem is well-posed. 



I. INTRODUCTION 



Current numerical relativity codes designed to simu- 
late the inspiral and merger of a black hole binary are 
limited by instabilities. These are now believed to be 
usually instabilities of the continuum equations, rather 
than of the numerical method. Intuitively one feels that 
the initial value problem for the Einstein equations is 
well-posed in a geometrical sense. For instance, in fully 
harmonic spacetime coordinates the Einstein equations 
reduce to ten quasilinear wave equations Q]. However, in 
the usual 3+1 split one uses only six of the ten Einstein 
equations for evolution, while the other four must be im- 
posed as constraints on the initial data. An evolution in 
which the constraints are not obeyed is not a solution of 
the Einstein equations. 

Consider perturbing a solution of the evolution equa- 
tions around a solution of the full Einstein equations. 
We call a linear perturbation physical if it obeys the con- 
straints initially and therefore at all times, and if it can- 
not be removed by a change of coordinates. A perturba- 
tion is called pure gauge if it obeys the constraint but can 
be removed by a coordinate transformation. A perturba- 
tion is called unphysical if it violates the constraints. All 
of these perturbations can be instabilities in the (weak) 
sense that they grow with respect to the background so- 
lution. An example for a physical instability is one that 
pushes a marginally stable star over the edge of gravita- 
tional collapse. An example of a gauge instability is given 
by the evolution of Minkowski spacetime with constant 
lapse, with an initial slice that is not quite flat: the time 
slices become singular in a finite time. Some, perhaps all, 
gauge instabilities can be avoided by a suitable choice of 
the lapse and shift adapted to the solution. 



The problem in numerical solutions are unphysical in- 
stabilities. These will always be triggered in a numerical 
simulation by hnite differencing or round-off error. Fur- 
thermore, the solution space of the evolution equations 
is infinitely bigger than the solution space of the full Ein- 
stein equations, and so the latter are likely to represent 
an unstable equilibrium. Evolution schemes which re- 
place between one and four of the evolution equations 
with constraint equations (constrained evolution) do not 
fundamentally address this problem, as they are still only 
solving six equations. 

In discussing instabilities, it is important to distin- 
guish between evolution systems which are well-posed 
and evolution systems which are ill-posed. In the for- 
mer, the growth of any linear perturbation Su(x, t) of a 
background solution uo(x,t) can be bounded as 
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where f(t) depends on uq but is independent of Su(x,0). 
This means that the solution uq + Su depends contin- 
uously on its initial data. By contrast, in an ill-posed 
system no such bound f(t) exists. Rather, the growth 
rate increases unboundedly with the highest spatial fre- 
quency present in 8u(x,0). These are instabilities in a 
rather stronger sense. They can be gauge, but typically 
are constraint- violating. 

A numerical simulation inherits all the instabilities 
already present in the continuum. A good numerical 
scheme does not add any others, but can never fix contin- 
uum instabilities. A crucial point is that in a numerical 
evolution, the highest spatial frequency present is effec- 
tively always the grid frequency, because it is generated 
by finite differencing error even if the initial data are 
smooth. 

As one increases the resolution in numerical solutions 
of an ill-posed system, constraint violating instabilities 
start with smaller amplitude (because they are initial- 
ized by finite differencing error), but grow more rapidly 
(because the finer grid can represent a higher spatial fre- 
quency). Therefore they don't converge away. At high 
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enough resolution and late enough time this will become 
apparent as a breakdown of convergence for the entire 
solution, but one can detect their presence before that 
by looking at a Fourier transform in space 3]. As these 
non-convergent instabilities are features of the continuum 
system they cannot be suppressed by any consistent nu- 
merical dissipation. 

By contrast, the numerical solution of a well-posed 
scheme will in general still be swamped by constraint- 
violating instabilities at late time, but now these start 
with smaller amplitude and grow at the same rate as one 
increases resolution. They therefore converge away. 

In practice, depending on the system, physical solu- 
tion, runtime and resolution, either convergent or non- 
convergent instabilities can be the dominant source of 
error. In 3-dimcnsional simulations in particular, the 
available resolution is quite limited, and the lack of 
convergence may not become apparent. The Arnowitt- 
Deser-Misner (ADM) formulation of the Einstein equa- 
tions has for many years been the main formulation used 
in 3-dimensional simulations, even though it is ill-posed 
(weakly hyperbolic) 2] . More recently a systematic com- 
parison of high-resolution, long-time evolutions for well- 
posed and ill-posed (weakly hyperbolic) systems Q has 
demonstrated that the breakdown of convergence is really 
inevitable, and has revived interest in well-posed formu- 
lations of the Einstein equations. 

This interest has focused on first-order reductions of 
the Einstein equations, because for general first order 
evolution systems useful criteria for well-posedness are 
known P|. A sufficient and necessary criterion for the 
initial value problem (with no boundaries, or periodic 
boundaries) to be well-posed is strong hyperbolicity: this 
roughly means that the system has a complete set of char- 
acteristic variables with real speeds. Symmetric hyper- 
bolicity is another criterion which implies strong hyper- 
bolicity and can be used to obtain a well-posed initial- 
boundary value problem: roughly speaking it means that 
the principal part of the system admits a conserved 
energy. Therefore a number of strongly or symmet- 
ric hyperbolic first-order reductions of the ADM equa- 
tions have been suggested over the years to assure well- 
posedness, see for example [H H 0, 0-In a symmetric 
hyperbolic formulation of the Einstein equations, the 
constraints can be dealt with consistently in the initial- 
boundary value problem (see for the full, and 0, Hl| 
for the linearised Einstein equations). 

However, first order reductions introduce new, auxil- 
iary, variables and constraints, and so further increase the 
solution space. One would expect this to give rise to addi- 
tional constraint- violating instabilities of the convergent 
type, even if well-posedness rules out the non-convergent 
type. There is some evidence that this is a real problem 
[12| , which may even outweigh the benefits of hyperbol- 
icity. It would therefore seem preferable to find a system 
that is symmetric hyperbolic while enlarging the solution 
space as little as possible, and in particular this should 
be a second-order system. 



In a companion paper |13| we have proposed a defi- 
nition of symmetric hyperbolicity for second-order sys- 
tems. This can be used to obtain a well-posed initial- 
boundary value problem in the same way as in first- 
order systems, but without enlarging the solution space. 
Here we show symmetric hyperbolicity for two formu- 
lations of the Einstein equations that are variants of 
the Arnowitt-Deser-Misner (ADM) equations, namely 
the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) sys- 
tem [Til Ejj, which is already widely used in numerical 
relativity El [H H fTgj . and a simpler related system 
suggested by Nagy, Ortiz and Reula (NOR) @- 

We know of two previous partial results concerning the 
hyperbolicity of a second-order, ADM- like version of the 
Einstein equations. In |20j| it was shown that a first-order 
reduction of the BSSN system is strongly hyperbolic, and 
a variant with some superluminal characteristic speeds is 
symmetric hyperbolic. It was then noted that the aux- 
iliary constraints associated with the introduction of the 
first-order auxiliary variables form a closed subsystem of 
the constraint system. This means that if the auxiliary 
constraints are obeyed initially, and if suitable boundary 
conditions are imposed, they are obeyed during the evo- 
lution, even if the other constraints are not obeyed. In 
this sense, the introduction of the auxiliary variables has 
not enlarged the solution system. In 0] it was shown that 
the NOR system (in second-order form), and its associ- 
ated constraint system are strongly hyperbolic, by using 
a pseudo-differential reduction to first order which also 
does not enlarge the solution space. This method relies 
in an essential way on Fourier transforms, and cannot be 
used to deal with the initial-boundary value problem. 

Here we go beyond these two papers in showing sym- 
metric hyperbolicity for the second order BSSN and NOR 
systems and their associated constraint systems, without 
any reduction to first order. We do this by defining char- 
acteristic variables, finding a conserved positive definite 
covariant energy, and expressing it in terms of the char- 
acteristic variables. We have chosen BSSN and NOR 
because BSSN is popular with numerical relativists to- 
day, whereas NOR seems a simpler version of BSSN that 
shares all its advantages. 

The paper is organised as follows. In Sec[H]we intro- 
duce our method and general notation. In Sec. II I II we 
prove symmetric hyperbolicity for the BSSN evolution 
equations and the associated constraint system, and in 
Sec. HVI we do the same for the NOR system. We take a 
first look at the initial boundary problem for both sys- 
tems in Sec. We propose a family of boundary condi- 
tions for a smooth boundary where the shift is tangential 
to the boundary. We show that there are no arbitrarily 
rapidly growing modes — this amounts to a strong indi- 
cation of well-posedness but is short of a proof. Sec. IVII 
contains our conclusions. 
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II. METHOD AND NOTATION 

Here we briefly summarise the relevant notation and 
methods of [l^. We need to split 3-tensors into their 
longitudinal and transversal parts with respect to a given 
direction Uj. Assume a 3- metric, say 7-y, which will be 
used to raise and lower indices, and let n l be a unit vector 
with respect to this metric. Then 

qj = - (2) 

is the projector into the space transversal to n.;. A tensor 
index that has been projected will be denoted by the 
index A, B, . . . instead of . . .. An index n denotes a 
tensor index contracted with ni or n l . A self-explanatory 
example of this notation is 

p i d l = p n d n + p A d A . (3) 

A pair qq of indices indicates a contraction with <jry , and 
for tensors of rank 2 or higher, we use the convention 
that they are totally tracefree on their projected indices. 

Now consider a system of evolution equations that are 
second order in space and first order in time. Linearise 
around a solution, and approximate the background- 
dependent coefficients of the linearised system as con- 
stant (frozen) . Retain only the principal part of the equa- 
tions. We call the resulting linear system with constant 
coefficients strongly hyperbolic if for any given direction 
n % we can find a complete set of characteristic variables 
U that obey 

dtU = Xd n U + transversal derivatives. (4) 

Here —A is the propagation speed in the rij direction. 
The U are constructed from u = (u, diu). This definition 
has two important consequences: transversal derivatives 
8au are automatically zero speed variables, and arbitrary 
multiples of transversal derivatives Dau can be added to 
any characteristic variable. Below we shall find the char- 
acteristic variables U in two steps: we first find a set 
of non-zero speed variables U' that do not contain any 
transversal derivatives, then add transversal derivatives 
to them until we can express the energy and flux in terms 
of the modified characteristic variables U. 

We call the system symmetric hyperbolic if it admits 
an energy 

E= I edV, (5) 
Jn 

where e is covariant, positive definite, and conserved in 
the sense that 

d t e = d t F l (6) 

for some flux F l , and if we can express e and F n in terms 
of characteristic variables. 

The simplest example is the wave equation in the form 

d t <j> = n, (7) 
d t U = ditffa (8) 



with u = (<fi, n) and u = (di4>, U). The characteristic 
variables in a given direction rii are 

u± = n±a„0, (9) 

U A = d A & (10) 

with speeds A = (±1,0). The covariant energy and flux 
are 

e = U 2 + d l <j>d l c/), (11) 
F i = 2n<9>. (12) 

In terms of characteristic variables they are 

e = \{Ul + U 2 _) + U A U A , (13) 
F n = \ {U 2 + -U 2 ). (14) 

Therefore 




where U± are now the characteristic variables normal to 
the boundary <9f2. 

The maximally dissipative boundary condition 

U+ = kU- + f (16) 

with \k\ < 1 and / a given function then guarantees that 
the growth of E is bounded by / and that E does not 
grow for / = 0. 

To make the energy positive definite in (f> itself rather 
than just di<f>, one can add a term a 2 (j> 2 to e, with a > 
constant. With a maximally dissipative boundary condi- 
tion and / = 0, E is then bounded by E(t) < e at E(0). 

Finally, allow the linearised system to have variable 
coefficients (from the non-constant background solution) 
and a non-principal part. For a finite time interval, the 
energy is then still bounded as E{t) < Ke at for some 
constants K and a, and the linearised initial-boundary 
value problem remains well-posed. Therefore it is suf- 
ficient to establish well-posedness to examine the prin- 
cipal part in the frozen coefficient approximation. On 
an intuitive level, this is so because the purpose of well- 
posedness is to rule out instabilities of the non-convergent 
type, which have a growth rate that grows with spatial 
frequency Well-posedness of the linearised problem is 
a necessary condition for well-posedness of the full non- 
linear (quasilinear) problem. 

III. THE BSSN FORMULATION OF THE 
EINSTEIN EQUATIONS 

A. Field equations 

The BSSN formulation of the Einstein equations (with- 
out matter) is obtained from the ADM form of the Ein- 
stein equations |2l| . 

dtjij = £/37ij -2aKij, (17) 
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d t Kij = CpKij - DiDja 

+a(R ij -2K il K l j +KK ij ), (18) 

H = R— A ,, A ' *' + AT 2 = 0, (19) 

Mi = D j K i i-D i K = 0, (20) 



by introducing the new variables 

lij = (det7)" 1/3 7y, 
1 = 7 J 7 7 ife j, 

= ^lndet7, 

N-l/3 



Ay = (det 7 )-^ \K ij -- lij K 



(21) 
(22) 

(23) 
(24) 



In the remainder of this Section, indices are moved with 
7y and its inverse 7 y . Generalising the BSSN equations, 
we densitise the lapse a with the determinant of the met- 
ric as 



(25) 



where a is a constant and now Q, rather than a, is a 
given function of the coordinates. The definition of the 
r l gives rise to the differential constraint 

G* = 7«r J - 7**7<i,fc = 0. (26) 

The definition of Ay gives rise to the algebraic constraint 

T = fiA l3 = 0, (27) 

and from the definition of 7^ we have the algebraic con- 
straint 



D = lndet7 = 0. 



(28) 



The BSSN equations are first order in time, second 
order in space, and quasilinear. The principal part of 
the evolution equations for Ay, K and T l is given by 
the highest spatial derivatives, (d 2 (j>, 9 2 7, dA, dK, dT). 
The evolution equations for 7 and <j) do not contain any 
of these highest derivatives. We define their principal 
part to be given by the next highest derivatives, that is 
(dcf>, c?7, A, AT, r). Note that with this definition of the 
principal part the equations are still quasilinear, because 
the evolution equation for 7^ is linear in K%j. Note also 
that f3 k dk is part of the principal part of the equations 
while d[3 terms are not. We define the derivative operator 



d = a-\d t - (3 k d k ). 



(29) 



It is the derivative along the unit vector field normal to 
the slices of constant time. The principal part of the 
evolution equations is then 



d Q (j> 



6 ' 
-2 .Aii' 



(30) 
(31) 



d K 
doA.j 



-6ere 

-40 



(32) 



I jij.mn 



2(1 + 3a)0,, 



TF 



r 



d r ~ 2(6 



l^^A^-^Kj, 



(33) 
(34) 



where ~ means equal up to non-principal terms, and TF 
indicates the tracefree part. We have added the con- 
straints aGnj), 2bMi and — (c/6)e -4 ^ 7^7^-0^ to the 
field equations with free coefficients a, b and c. Adding 
these terms changes the evolution off the constraint sur- 
face which can affect the hyperbolicity of the system. 

The principal part of the Hamiltonian and momentum 
constraints is 



H = H 

Mi 



a Gj,j ^ 



0. 



(35) 
(36) 
(37) 



Here a' and d parameterise different ways of writing the 
Hamiltonian constraint that are found in the literature. 
We shall work explicitly only with Hq, and so a' and c' 
will not appear below. There are many versions of the 
BSSN equations which vary in small details in both the 
principal and non-principal parts. For comparison, the 
principal part of the version given in [2^ is characterised 
by ct = (the lapse is not densitised) and a = b = a' = 
c' = l, c = 0. 

The constraints are compatible with the evolution 
equations, which means that they form a closed evolu- 
tion system. It is 



d H 
d Mi 

doG, 
d T 
d D 



-2e-^f j Mi d , 
\H , i + e-^f k Gi, jk 



V2 

+-^ k G j<ik + — 
2bMi, 



V k D 



-T 



i]k 



2 

-2T. 



(38) 

(39) 
(40) 
(41) 
(42) 



B. Strong hyperbolicity of the main system 

For the purpose of decomposing 3-tensors and tensor 
equations, we define m, n l and qij with respect to the 
conformal metric 7™. In the frozen coefficients approxi- 
mation, the undifferentiated 3-metric (7^ , j 13 and <fi) is 
to be treated as a background quantity, while Ay, T l , 
K, 7ij fc, 4>,i and their derivatives are to be treated as 
dynamical variables, and decomposed with respect to n,. 



5 



We now prove strong hyperbolicity of the main and 
constraint systems by constructing a complete set of char- 
acteristic variables U. The quantities 



u = (di4>,di^ jk ,A i: j,f\K). 
obey the pseudo-first order system 

d t u ~ (aP z + F)diU 

or equivalently 

d$u ~ P l diU. 



(43) 



(44) 



(45) 



(This is not a genuine first-order system because in ex- 
pressions like di(dj(f) that appear on its right-hand side 
we allow ourselves to commute partial derivatives, rather 
than treating dj<f> as a 1-form variable dj, as we would in 
a genuine reduction to first order.) As discussed in 
transversal derivatives are automatically zero speed char- 



acteristic variables. Here we have P n dA<t> = Pn^Ali 



0. 



We shall write the eigenvalues of P n as Ae -2 *. With this 
definition A measures the propagation speed in units of 
the speed of light, measured with respect to the do ob- 
servers. In particular A = ±1 variables propagate along 
the light cone. 

In tensor components with respect to rii, the scalar 
block of P n is given by 



(46) 

(47) 
(48) 
(49) 



P n d n <f> = 


- l -K, 

6 ' 




—94 






PnK = 




PA - 
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1 — C 



(50) 



P A 

1 n-fiqq 



3 ^ n ^n^fnn) 



6 



l + 2c 



G 



p r — 
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+ 2(6-1) A nn . 



(51) 
(52) 



The two algebraic constraints T = and D — are often 
enforced in numerical simulations after each time step. 
We can mimic enforcing T = in our analysis by setting 
A nn +Aqq = and dropping one of the two variables from 
the system. Similarly, we can mimic enforcing D = by 
setting <9„7 nrl -I- <9„7 9g = and dropping one of the two 
variables from the system. Both reductions affect the 
hyperbolicity of the system. 



a. Algebraic constraints not enforced The scalar 
block is diagonalisable for 



c > 0, 77 > 0, (7 > 
with eigenvalues 

where we have defined the shorthand 

4a6- 1 
V = — o — ■ 



(53) 
(54) 

(55) 



(Allowing the system to be strongly hyperbolic without 
enforcing the algebraic constraints is the reason why we 
have introduced the c term, which potentially gives doAij 
a non-zero trace.) 

b. Trace constraint enforced If we enforce T = but 
not D = this is consistent with the evolution equations 
only for c = 0. The scalar block is diagonalisable for 



with 



0, 77 > 0, a > 0, 



A= {0,0, ±^77, ±v^} 



(56) 



(57) 



It would be inconsistent to enforce D — but not T = 0. 

c. Both algebraic constraints enforced If we enforce 
both algebraic constraints D = and T = the scalar 
block is diagonalisable for 



rj > 0, cr > 0. 



(58) 



(Note that c becomes irrelevant in this case.) The char- 
acteristic variables that do not contain any transversal 
derivatives are: 



U' Q = 
Ui = 



with speeds 



(b - l)d n % n + f n - 8bd n (j), 
(1 - 4a)<9„7„„ + 4af„ - 8d n cj) 
±e 2 *JT 1 {&A nn -AK), 



The vector block of P n is 



Pn dnlAn 
PtlAati 
Pnf A 



-2 A 

-40 



(f A - d r 



7-Ar 



2(6-1)4, 



(59) 

(60) 
(61) 

(62) 

(63) 
(64) 

(65) 



as well as the trivial P n dA<P — P n dAlnn — PndA%q = 0. 
It is diagonalisable for ab > 0, with nontrivial character- 
istic variables 



U' A 

U' ±A 



(b-l)d n % A + f A , (66) 
-a<9«X4„ + of a ± 2\/^e 2 ^i A „ , (67) 
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with speeds A = (0, ±Vab). The tensor block is 



The vector sector of the constraint propagation is 



PndnlAB = ~2A AB , 

P n A AB = -e~ 40 i<9„X4B, 



(68) 
(69) 



as well as the trivial P n d A jBn = 0. It is always diago- 
nalisable, with nontrivial characteristic variables 



(70) 



with speeds A = ±1. We see that the vector and tensor 
sectors do not add any new conditions for strong hyper- 
bolicity. 

In summary we have shown that the BSSN system 
with both algebraic constraints enforced continuously is 
strongly hyperbolic if a > and rj > 0. If neither al- 
gebraic constraint is enforced, c > and c ^ r/ are also 
required. On the other hand, if only the trace constraint 
is enforced, c = is required. (If both constraints are 
enforced, c is irrelevant and can be set to zero.) The 
characteristic spectrum of the complete system is 



A = {0, ±y/a, ±y/rj, ±Vab, ±1}. 



(71) 



C. Strong hyperbolicity of the constraint system 

We construct characteristic variables for the constraint 
system from the set 

u c = {H , M u diG^diT, didjD). (72) 

The scalar sector of the constraint propagation is 

P n Ho = -2e-^M n , (73) 

P n M n = ±H + e~ 4 * (y 0»G„ + (74) 

P n d n G n = 2bM n , (75) 

P n d n T = - C -e~^dlD, (76) 

P n d 2 n D = -2d n T. (77) 
If we dehnc 

C' ± = Aad n G n ± 6V?ye 20 M„ + e 4 ^H (78) 

the characteristic variables are 

C Q = e^bHo + d n G n , (79) 
\/cd 2 D =f 2e 2 ^d n T, (80) 

(1 - c) [d 2 n D T ^ e209 « T ) + (l - £) C '±> ( 81 ) 

with speeds A = (0, iy/c, iy/rj). If we optionally impose 
the trace constraint, the characteristic variables are C 
and C'± + d 2 D with speeds A = (0, ±y/rj)- If we impose 
both the trace and determinant constraints, we are left 
with Cg and C' ± with speeds A = (0, ±y/vi)- 



P n M A = e - 4 ^d n G A , 
P n d n G A = 2bM A , 



(82) 
(83) 



as well as P n d A G n = P n d A T = P n d A d n D = 0. The 
nontrivial characteristic variables are 



CL A = ad n G A ± 2y/abe 24 'M A 



(84) 



with speeds A = ±y/ab. The tensor sector is completely 
trivial, with P n d A GB = P n d A dBT = (including the 
traces). The constraint system, in all three cases, is 
strongly hyperbolic as long as the main system is strongly 
hyperbolic. 



D. Main system energy 

We look for an energy density e of the second-order 
system that is positive definite in Aij, K, l\, 7^ and 
(j>.i and obeys a conservation law. With the shorthands 
U = J jk Jjk,i = D,i and d, 
quadratic form in these variables is 



7 jfe 7ii,fe, the most general 



e = c e 4 ^ A i:j 2 + a Ti 2 + c 2 t 2 + c 3 di 2 

+c 4 djf * + c 5 Ur + c 6 dit* + c 7 e 4 ^ K 2 

+c 8 7y,fe 2 + c 9 7u,fe7 ifeJ + cio f ^ + cue 40 TK 

+ci 2 4>/ + C13 d^ 1 + c u UV + cise 40 T 2 . (85) 

In the cross term jij i k'J lk ' : ' and similar terms in the re- 
mainder of the paper indices are raised only after differ- 
entiation. The general ansatz for the flux F l contains 14 
free coefficients. After some linear algebra it is possible 
to see that this e and F l together obey a conservation 
law only if 



4ab = 1 + 3(7, 



(86) 



or equivalently r\ = a. 

There are additional restrictions depending on whether 
or not we enforce the algebraic constraints during evolu- 
tion. If we work with non-vanishing D and T then we 
need ab = a = 1 and cn(c — 1) = 0. On the other hand, 
if we enforce both algebraic constraints the parameter c 
and several of the coefficients Ck become irrelevant, with 
no additional constraints on the parameters. 

From now on, in order to simplify the equations, we 
enforce both the T = and D = constraints in the 
remainder of this Section. 

The generic conserved energy depends on four free co- 
efficients, one of which is an overall factor. With the 
shorthands 



= e 4 W + + 



2a — ab — 1 _ 



-7ij,k7 



-adi(r 



(87) 
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£2 
£3 



r, - 



+ (&-l)dj 



e 4 ^ 2 + 36<r0' ; 



(88) 

(89) 
(90) 



the energy density is 

e = c e + cjei + [c 3 - (6 - l) 2 ci]e 2 + c 7 £3- (91) 

The free parameters Co, ci, C3 and C7 will be restricted 
below by inequalities derived from the requirement that 
e be positive definite. 

This energy reduces to that given in J2jj for a first- 
order reduction of BSSN with the change of notation 
a —> 2a, b —> m, % jik -> d kij and 4> ti -> d.,/12, 
the restriction a = 1 to the parameters of the evolu- 
tion equations, fixing the overall scale as Co = 1, and 
the further restrictions c\ — 1/(26 — 2), C3 = and 
C7 = l/a = 3/(46 — 1) on the coefficients of the energy. 
The condition l|86[) on a reduces to a similar condition 
in [2(j. The choice of Sarbach et al. can be interpreted 
as follows: C3 = and c\ = a/ (2b — 2) eliminate di from 
the energy (03 = c± = C13 = 0); assuming those condi- 
tions, a = 1 is the only choice that eliminates the cross 
term 7ij j fc7 fcl '-'. Their choice for c 7 has no relevant effect. 
Note that the positivity condition on their choice of c\ 
forces b > 1, and together with a = 1 this gives rise to a 
super luminal speed. By contrast, by retaining the con- 
tractions ti and di in the ansatz for the energy we shall 
be able to make all speeds physical. 

The four terms in the energy are conserved separately: 



dot! = 0, 
8 e 2 = d^l 
d e 3 = diFi 



(92) 
(93) 
(94) 
(95) 



(Condition l|86|) is required only for the conservation of 
eo-) The fluxes are given by 

2A VJ [atj - a{b - l)dj - 2(1 + 3a)(j)j] 

+f l A jk [2(l + ab- 2<z) 7Afe - j jk ,i], (96) 

4A^ k (f m j mjtk - f 3 % m ' m ), (97) 

-12aK(j)'\ (98) 



Fi 



The energy density in terms of characteristic variables is 



e = %(Ul AB + U*_ AB ) + £-(Ul A + U 2 _ A ) 



c c 7 



Aab 

(ul + u 2 _) + ^±^(v 2 



16cr(2co + 3c7)^ + 6 
+ quadratic in zero speed variables 

and the flux is 

T 2 tt-2 ^ , C 



■V 2 ) 
(99) 



ftp" = ^ {u : +ab 



- u- 



AB 



) + ^h(U 2 +A 



C C 7 



16^(200 + 3c 7 ) 
2c + 3c7 i— fj 2 



4Va6 

(ul - U 2 _) 



u 2 . A ) 



6 



a (VI -VI 



(100) 



We have modified the characteristic variables by adding 
terms in d A "fij in order to write e in terms of character- 
istic variables. With the shorthands 



Z\ — 


c 3 -(6-l) 2 ci 








CO 


Z2 = 


1 + 3a6 — 4a + 6z±, 


(102) 


Z3 = 


a — ab — 2z\, 


(103) 


Zi = 


2 + 3a6 — 5a + 621, 


(104) 


Z5 = 


1 + ab - 2a + 2z\, 


(105) 


Z6 = 


CO 


(106) 


2^(2^ + 3c 7 ) ' 



the modified characteristic variables are 

U± = (1 - 4a)d n g nn + 4af „ - 8d n (f) 

±V^e 24, (6A nn - iK) - 2z 2 d A g An , (107) 

V± = Qy/ad^^e^K - z 6 U±, (108) 
U± A = ±2\fabe 2 ^ A An + af A - ad n g An 

+z 3 d B gAB + -^d A g nn - 8abd A <j>, (109) 

U± AB = T^Aab + 7}d n g AB 

-y (cUpsn + d B g An - q AB d c g C n)-(110) 

The expressions for Uq and U A are not given, but they 
also include transversal derivatives. Note that because 
rj = a, U± and V± are eigenvectors for the same eigen- 
values, and we have used this to make the energy diagonal 
in U± and V±. 



E. Constraint system energy 

The constraint energy is quadratic in Hq, Mi and Gij. 
(We must use Gi.j rather than Gi so that all terms are 
of the same order in derivatives). The most general con- 
served expression is 



woe^H 2 + \2w x e^M 2 



3awi 



+w 2 G iij CP> i + 2 ^±^e^G/H 



wq + (1 + ab)w\ — b 2 W2 
b 2 



{GSf 



(111) 



with arbitrary coefficients Wq, W\ and w<i . The corre- 
sponding flux is 

Fl = 4w 1 (e 4 ' t 'H M i + WaG-'M, + aAPG j d ) 

+4bw 2 (G l ' j Mj - U' G'',,). (112) 

In terms of characteristic variables, the constraint en- 
ergy and flux are 

e c = ^ C -AC A + C +A Ci) + ^{C 2 _ +C 2 ) 
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+ quadratic in zero speed variables, (113) 



e 2<t>pn 



3wi 



2Vab 



(C +A C* - C- A C*) 



(114) 



where the modified non-zero speed characteristic vari- 
ables are 



C± = e 4 ' t >H ±6e 2 *VvM n + Aad n a n 



(115) 



C±A 



w 2 b , 



±2v / o6e 2 *M A + ad n G A + —d A G n .(lW) 

3wi 



Note that w 2 does not appear explicitly in or (|114fl , 
but it does appear in the definition of the characteristic 
variables. 



F. Symmetric hyperbolicity and causal speeds 

With the condition r\ — a that is required for conser- 
vation of the main energy, the complete spectrum can be 
written entirely in terms of a: 



A 



VI + 3ct 



±1 



(117) 



Strong hyperbolicity requires a > and the absence of 
superluminal speeds requires a < 1. Note that with a = 
ab = 1 all speeds are either or ±1. 

Now we impose positive definiteness of the energies, 
that is, all eigenvalues of their respective matrices are 
strictly positive. This means that the main energy van- 
ishes if and only if 7^ and <j> are constant and all other 
variables vanish, and that the constraint energy vanishes 
if and only if Gi is constant and the Hamiltonian and mo- 
mentum constraints vanish. (We are assuming that the 
two algebraic constraints vanish because they are being 
continuously re-imposed.) 

To see when the constraint energy is positive definite, 
we need to decompose Gij as 



(118) 



where Sij is symmetric and tracefree and is anti- 
symmetric. We obtain a sum of simple squares and a 
quadratic form in the two variables Hq and G % 

The constraint energy is positive definite if and only if 



wq > 0, < wi < 3wq, —7— — - < ab < 1, 



-ab < 



w 2 b 2 
3wi 



< ab 



w + W\ 
2w 



(119) 



It is interesting to see how w\ affects the possible range 
for ab. 



The positivity of the main energy is more complicated. 
We must take into account that the partial traces ti and 
di of the three-index object 7^ appear both explicitly 
and explicitly in the energy. We therefore decompose it 
as 

lij.k = ■= [(3d(i - %)7j)Jt + ( 2i fc - dk)%] + fijk (120) 



where fijk is completely tracefree. We then obtain a 
quadratic form in ti and di, plus csffj k + cgfijkf lk -' ■ To 
analyse the positivity of the latter we decompose fijk 
in its 12 independent frame components and analyse the 
corresponding quadratic form by brute force. It turns out 
to be positive if and only if cs > and — cs < eg < 2cg. 
Finally, assuming causal speeds, we obtain the positivity 
conditions 

c > 0, c 7 > 0, 2a 2 c < (4a6 - l)ci, 
2a 2 c - (4afe - l)ci 



5c - 



36ci 



< C3 - (6- l) 2 ci + 



1 + 2ab - 4a 



c < 0. (121) 



The three sets of inequalities are compatible, and can 
be solved sequentially, in this order: first choose cq, C7, 
wo and w\ independently, then a and 6, followed by C\. 
Finally choose C3 and w 2 - It is clear from this construc- 
tion that the solutions form a single connected set. For 
example, a possible solution with A = (0, ±1) is 

a = 6 = cr = l, (122) 
00=07 = 1, ci=2, c 3 = 0, (123) 

u>o = l, w x = ~, w 2 =0. (124) 



IV. THE NOR FORMULATION OF THE 
EINSTEIN EQUATIONS 

A. Field equations 

In this Section, all indices are moved with 7^ and its 
inverse 7 y . The NOR system is obtained from the ADM 
system with densitised lapse by introducing the variables 



fi = lif 3 



P ik 

-jl Tjfc.ii 



(125) 



where p is a constant parameter. In 0] the choice p = 1 
is made, p = 2/3 makes fi — jijT J . The definition of fi 
gives rise to the constraint 



Gi = fi- jij' J + 



P jk 



0. 



(126) 



As in the BSSN system, we parameterise the use of this 
constraint with a and the use of the momentum con- 
straint with b by adding aGufl and 2bM; to the evolu- 
tion equations for Kij and fi. Following 0], we also add 
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cjijH to the evolution equation for Kij with free param- 
eter c. The principal part of the evolution equations is 



dolij 
doKij 



do f x 



-2K ijt (127) 

+ (1 - a){j klJ ' k +"f k jX k ) 

+ (ap - 1 - <7)7 H 7M,tfl + nijH, (128) 



2(6 -l)i^ 



(P - 2b)K t1 



(129) 



and the tensor sector is always diagonalisable, with A = 
±1. For general values of (p, a,a,b, c, a') the characteris- 
tic variables are too long to give here. 

The scalar sector of the constraint system is diagonal- 
isable for x > with characteristic variables 

C = 8 n G n + bH , (141) 
C ± = (1 + 4c)H + 4ca'd n G n T 2^M n , (142) 

with speeds A = {0, ±^/x}, and the vector sector is diag- 
onalisable for ab > with 



where 



H = Ho + a'Gi 



Ho 

Mi 



K ' j 



K A . 



In @, p = o = a' = 1. 
The constraints system is 



doHo 
d Q M z 

dad 



-2Ms\ 



1 



^2 

a ^ ,, 
2'" 
2bM,. 



2c H , 



^ + 2ca' ) G 



(130) 
(131) 
(132) 



(133) 

(134) 
(135) 



B. Strong hyperbolicity 

This proceeds exactly as in the BSSN system, with 
pseudo-first order variables 



and 



u= (dijj k ,Kij,fi) 



= {H ,M i) d i G j } 



(136) 



(137) 



for the constraint system. We define n.j, n % and qij with 
respect to 7^ . (In |2| , a flat auxiliary metric is introduced 
instead.) 

The scalar sector of the main system is diagonalisable 
for x > 0, <J > and x^ff, with 

A = {0,±VX,±^}, (138) 

with the shorthand 

X=l + 4c(l-a'&). (139) 

If 77 = x we can diagonalise the scalar sector also for 
X = o. (Note that while we use the same shorthand r\ 
defined by f53|> as in the BSSN system, ^/rj is not a speed 
in the NOR system.) 

The vector sector is diagonalisable for ab > with 



A = {0,±Vab}, 



(140) 



C±a = ad n GA ± 2v abMA, 



(143) 



with speeds A = {±vab}. 

The union of conditions for both the main system and 
the constraint system to be strongly hyperbolic, and for 
all speeds to be physical (|A| < 1), is 



< a < 1, < ab < 1, < x < 1, 



(144) 



and either x 7^ or X — a — V- The union of all speeds 
is 



A = {0, ±1, ±Va~b, ±V^, ±y/x}- 



(145) 



We see that we can make all of these speeds either zero 
or one by choosing 



ab = a = x = 1 => c(l - a'b) = 0, 



(146) 



that is, either c — or a' = a. From now on we restrict 
ourselves to the choice c = 0. This makes x = 1, and 
removes a' from the system. We keep a, 6, p and a free, 
with the single restriction that a = 1 requires ab = 1. 
The value of p is irrelevant for strong hyperbolicity. Our 
results on strong hyperbolicity of the NOR system con- 
firm and generalise those of |2|, which were obtained us- 
ing a pseudo-differential reduction to first order. 



C. Main system energy 

The most general ansatz for the energy density e is 
quadratic in "fij,k, and /j. As well as contracting the 
free indices on pairs of these, we can form the contrac- 
tions K, di = jij'i and U = ^ k ^ijk,i first: 

e = c Kij 2 + ci f 2 + c 2 t 2 + c 3 d 2 



+c 4 dif + c 5 UP + c & dit 1 + c 7 K 2 

+C8 Ji 3 ,k 2 + Cg 7y,fc7 lfe J • 



(147) 



The most general form that is conserved depends on four 
free parameters Co, ci, C3 and C7 obeying 



2(ab - l)c 7 = (a - 2ab + l)c . 



(148) 



for arbitrary (p,<j,a,b). The coefficient c must be 
strictly positive and therefore there are two possibilities: 
either ab — 1 (which now implies a = 1) and then Cq 
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and c 7 are independent, or ab ^ 1 and then c 7 is deter- 
mined by Co and the parameters (a, a, b). What follows is 
valid in both cases, with the corresponding restrictions. 
Note that strong hyperbolicity and energy conservation 
together give the two directions of u = 1 ab = 1 . 
With the shorthands 



U) 

£2 
£3 



9 1 o 2a — ab — 1 
*V + ^fc 2 + g 



7ij,fe7 



-ad,- 



>-('-!) 



/, - (fc - § J ti + (6 - 1)* 
2 + 2ab~ 2ap + a 2 



(149) 

(150) 
(151) 

+ (a-l)d i t i -aUf (152) 

(the indices on jij.k are raised only after differentiation), 
the energy can be written as 

e = c e + c\e\ + [c 3 - (6 - l) 2 ci]e 2 + c 7 e 3 . (153) 

The flux is 

co^ + Ics-^-lfd^ + cr^, (154) 



F' 



with 



F' 

1 n 



*3 - 
Fi = 



(155) 
(156) 



2aK i i[f j -(b-l)d j -(b-p/2)t j ] 
+Y l K* k [2(1 + ab - 2a) 7 y, fc - 7,-*,, 

*y a K* k (>y l j,k-'njdk), 

K [2a f + (ap - 2 - cr)f + 2(1 - a)cf ] 
+2(1 - ab)K iJ tj. (157) 

(As in BSSN, e-y has no flux.) The flux in terms of char- 
acteristic variables is 



pn _ C (tt2 

b - y( u +ab 
c + c 7 



- c/ 



CO 



-AB) 



(U 2 



2^ 



(V+ V -) + Tk 



Wab 

c c + 3c 7 



16 cq + c 7 



-A U_ A 



(158) 



The characteristic variables including transversal deriva- 
tives are 



U± 



-dr^lqq ~ %Z 2 d 7 An ± 2-ffgg, 

__ ap — a — 1 

a Jn V <J K nn H O n "fqq 

ap — a — 2a 



-Onjnn + Z 3 d "/An 



(159) 



(160) 



c? 



-(-d n j qq - 2z 2 d A jAn ± 2y/aK qq ), 



U± A 
U±AB 



2(c + c 7 ) 
o/a - ad n "fAn ± 2VaWdn 

+2; 3 5 B 7AB + dA(z±"/nn + Z 5 "f qq ), 
TKab + \<>nlAB 

+Z2 (d(AlB)n - \<lABd C ~1Cn 



(161) 



(162) 



with speeds A = {±1, ±\/cr, ±va6, ±1}, where we have 
defined the shorthands 



C3 _(6_l)2 Cl 



C() 



(163) 

(164) 
(165) 



Z2 = 2a — ab—l — 2z\, 

Z3 = a — ab — 2z±, 

Zi = (2a6-o--l)/2 + l-2a + ap/2 + 2zi,(166) 

z 5 = (2o&-a--l)/2 + a(l-3& + p)/2-«i.(167) 



Note that the last term of (|160|l requires Co + c 7 7^ 0. We 
shall see below that this is always true. 



D. Constraint system energy 

The constraint energy is quadratic in Hq, Mi and Gij. 
The most general form that is conserved has three free 
parameters wq, w\ and wi for arbitrary (p,a,a,b). It is 

e c = WoHq + AwiMf + ^p-(Gi,j) 2 

+w 2 G id G>> i + 2 Wo ~ Wl G\ i H Q 



w — (1 + a6)wi — 6 2 W2 
b 2 



{G\i) 2 (168) 



The flux is 

Fl = -4(awi + bw 2 )G\ ] M i - 4w 1 H M i 

+4{bw 2 G ij +aw 1 G j ' i )M j . (169) 

In terms of characteristic variables, 

Fl 1 = ^{C\ -Cl) + ^ b (C 2 +A ~ C 2 _ A \ (170) 
where, including transversal derivatives, 

C± = H +(a+ ^ d A G A T 2M„, (171) 
C ±A = ad n G A + h —d A Gn±2^/al>MA, (172) 

Wl 



with speeds A = {±l,±Vo&}. 



E. Symmetric hyperbolicity and causal speeds 

Positivity of the constraint energy requires that for a 
given value of ab we must choose Wq, wi and w 2 obeying 



wi 2 w 2 _ 3 l ^ wi 



w > 0. < — < 1, b z — < - 1 

w wi 4 \ wo 



^ <a6 <r,i 



(173) 

2 V w J wi 



The positivity conditions of the main energy will be anal- 
ysed separately for the generic and special cases. 
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In the general case ab ^ 1 and a ^ 1 positivity of the 
main energy requires that (using the strong hyperbolicity 
condition a > and the causal speeds condition a < 1) 



c > 0, 



1 + 3(7 



< ab < 1, 



a 2 Co(l — c) < 2erci(l — a&), 



5c 



2crci(l — afr) — a 2 c (l — cr) 



4 2ci(l - a6)(4 - Aab + 3<r) - a 2 c (4a& - 1 - 3cr) 

(174) 



l + 2afe-4a 
< c 3 - (6 - l) 2 ci H c < 



The value of p is irrelevant also for symmetric hyperbol- 
icity. The restrictions on a and ab guarantee that Co + c 7 
and Co + 3c 7 are always strictly positive. A simple exam- 
ple is 



a = 
co 



3 

T 



6 = 1, 



Cl 



w = 1, twi = 



1 2 

CT= 2' P= 3' 
c 3 = c 7 = 0, 

-, W2 = 0. 



(175) 
(176) 
(177) 



In the special case ab = a = 1 positivity of the main 
energy requires 

Co > 0, c + 3c 7 > 0, b 2 c\ > c + c 7 > 0, 
5cq b 2 ci — Cq — c 7 



4 36 2 ci — Co — 3c 7 
< c 3 - (6 - l)^ci H t^cq < 0. 



A simple example is 

a = b = 
co = 1, 
w 



Ab 



ci = 3, c 3 = c 7 = 0, 



1, wi = -, w 2 = 0. 



(178) 



(179) 
(180) 
(181) 



V. CONSTRAINT-PRESERVING BOUNDARY 
CONDITIONS 

A. The boundary system 

We have proved symmetric hyperbolicity for both the 
main system and the constraint system of the BSSN and 
NOR formulations. Therefore both formulations admit 
constraint-preserving boundary conditions of the type 
proposed in [ljj. A general discussion is left to future 
work. Here we summarize the basic idea, and propose a 
family of boundary conditions for the case of a smooth 
boundary with tangential shift. 

Because the constraint system is compatible with the 
evolution system, for every pair of characteristic variables 
C± of the constraint system with speeds ±A, there is a 



pair of characteristic variables U± with the same speeds, 
such that (after suitable normalisation), they obey 



C± = d n U ± 



(182) 



where the dots, here and in the remainder of this sub- 
section, indicate transversal derivatives and lower-order 
terms. In order to guarantee that the constraint en- 
ergy does not grow, we formally impose the homogeneous 
maximally dissipative boundary condition 



C+ - kC- = 



(183) 



on the constraint system. This must then be translated 
into a boundary condition on the main system. From 
(J2U, we have 



cW± = (±\a + f3 n )d n U± 



(184) 



In the following we restrict consideration to the case 
where (3 n — on the boundary, so that the A = char- 
acteristic variables propagate along the boundary. I|183|) 
is then equivalent to 



d t U+ + nd t U- = 



(185) 



We define a variable X = U+ + kU- that is restricted to 
the boundary, with evolution equation 



d t X = . . . , 
and impose the boundary condition 

U+ + kU- = X 



(186) 



(187) 



on the main system. 

We have chosen our notation in the previous sections 
so that in both BSSN and NOR we have 

C± = d n U± + ..., (188) 
C ±A = d n U±A + ..., (189) 

and the boundary conditions for either system are 



C+ - kiC_ = 0, 

C+A — K2C-A = 0, 
V+-k 3 V- = F, 
U+ab — K4U-AB — Fab, 



(190) 
(191) 
(192) 
(193) 



where F and Fab are free boundary data, and I|190I191[) 

are implemented as 



U + + kU- = X, 
U+a + hU-a = X A - 



(194) 
(195) 



In general the boundary system is coupled to the bulk 
system, so that X and Xa in (|194I195|I cannot be consid- 
ered as given a priory, and then the constraint-preserving 
boundary conditions are not true maximally dissipative 
boundary conditions. There are two exceptions, which 
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have been called "Neumann" and "Dirichlet" boundary 
conditions in the literature, where an extended bound- 
ary system can be given that decouples from the bulk 
system. This happens for the Einstein-Christoffel sys- 
tem linearised around Minkowski spacetime, |10| . the 
full Einstein equations in harmonic gauge |ll| . and the 
Maxwell equations fl3| . This boundary system can then 
be evolved before the bulk system is evolved, X and 
Xa can be treated as given a priori, and the constraint- 
preserving boundary conditions become true maximally 
dissipative boundary conditions. The details, assuming 
a smooth boundary with tangential shift, are given in 
Appendices 151 and IO However, we would like to stress 
that these boundary conditions are very restrictive: it 
does not seem very physical to find the boundary of the 
spacetime without knowing what is inside. 

As we have assumed that the shift is everywhere tan- 
gential to the boundary, and this is possible in the case 
of a non-smooth boundary (for example a cube) only for 
zero shift, we also restrict the analysis in this paper to a 
smooth boundary. 



B. Mode analysis 

If the boundary system does not decouple, we cannot 
use our current energy estimates to prove well-posedness 
of the initial-boundary value problem. We can, however, 
check a necessary condition for well-posedness, namely 
that there are no modes that grow exponentially in time 
where the growth rates increases unboundedly with spa- 
tial frequency. We conjecture that this condition is also 
sufficient. 

In the frozen coefficients approximation that we have 
been using throughout this paper we assume that the 
linearised perturbation varies over space and time scales 
much smaller than those given by the background so- 
lution and the numerical domain. For consistency we 
must therefore assume that the domain is a half space 
and that the boundary is a plane. We introduce co- 
ordinates so that the domains is — oo < x 1 < and 
— oo < x A = (x 2 ,x 3 ) < oo, and the metric in these co- 
ordinates is Sij. In the frozen coefficient approximation 
a > and j3 l are also constant in space and time. (As 
before we assume /3„ = on the boundary, and therefore 
everywhere in the frozen coefficient approximation). 

After a Fourier transform in x A and Laplace trans- 
form in t we are left with a system of linear ODEs with 
constant coefficients in x . In general this can be trans- 
formed into a matrix eigenproblem by an exponential 
ansatz in x 1 . (A priori an exponential times a polynomial 
in x 1 could be required to find the general solution but 
this is not the case here.) The general solution with ho- 
mogeneous boundary conditions can therefore be written 
as a sum of modes of the form 



Dqu = su, d n u = fiu and Oau — iloau. 

If the initial data are nonzero, we must add to the sum 
over modes a particular integral that does not concern 
us here because it is controlled by the boundary data 
Q. We are interested only in modes that are square- 
integrable over space at any moment in time. Therefore 
we assume that lua is real, and that Re [i > 0. s and fi will 
in general be complex. If a mode of this form exists for 
some (s, u>a, n, u), then one exists also for (ks, ku>A, k/J,, u) 
for any k > 0. Therefore, if any growing mode, with 
Re s > 0, exists, there are growing modes with arbitrarily 
large growth rates and the problem is ill-posed. A neces- 
sary condition for well-posedness of the initial-boundary 
value problem is therefore that the homogeneous bound- 
ary conditions rule out the existence of any mode with 
Re s > for real u>a and Re fj, > 0. 

For simplicity we concentrate again on NOR with (|179l - 
11811) . u — (jij, Kij, fi) is decomposed with respect to 
the normal vector = (1,0,0). It is helpful to intro- 
duce the notation f m = iu A fA and f p = p l fi where p 1 
is orthonormal to uji and rii, and similarly for other ten- 
sor components. Substituting the ansatz l|196[l into the 
NOR evolution equations, we find after some linear alge- 
bra that, for s ^ 0, 



(V 



and 



Ki. 



- M 2 )7y = 



fp = o, 

2 

fm = -gW 2 (7„„ +% q ), 
2 

fn = ^{inn+lqq)- 

For a non-zero solution to exist, we must have 

2 2 2 

H — sr = lo . 



(197) 

(198) 
(199) 
(200) 

(201) 
(202) 



u(x 1 , t) 



e (as+if3 A LUA)t+iuJAX A +iJ,x 1 i 



(196) 



where u is a constant vector. With this ansatz we have 



The coefficients %j are then free parameters. They de- 
termine the coefficients Kij and fi through (|198I201|I . 

Similarly, for the BSSN system with (|122I123|I we find 
the equivalent of condition (|197|l for 7^ and <j>, and 

(203) 
(204) 

(205) 
(206) 
(207) 

Note that sj^ represents the time derivative of 7^. 
Therefore, as we assume that the algebraic constraint 
D = is being imposed continuously, 7y must be trace- 
free. Similarly, from the algebraic condition T = 0, Aij 
is also tracefree. With (|197fl obeyed, 7^ and <j> are free 

coefficients, and determine Ay, K and T l . 



K 


= —6s4>, 




s - 








= 0, 


f 


= -8u; 2 4> 


fn 


= 8^0. 
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We now substitute the ansatz l|196(l with these coef- 
ficients into the six constraint equations i|19UH93l) . We 
obtain six algebraic equations that are linear in the four 
components % n , % m , % p , % q of % and the two com- 
ponents, 7 mm = — 7 P p, 7m P of the tracefree transversal 
object jab- We can solve these recursively to find that 
all 7^ = 0, as long as 

(1 - Ki)n + (1 + Ki)s ^ 0. (208) 

for all four Ki. For a mode to exist, this inequality must 
be violated for at least one of the Ki- Let the value of 
this Ki be k, which therefore obeys 

K (fj,-s) = (fi + s). (209) 

We now investigate the space of possible solutions 
([i,s,uj,k) of the two algebraic equations l|202|l and l|209|) 
with Re /i > and Re s > 0, with the aim of finding a 
condition on k that excludes all such solutions. We first 
consider the case u> = 0. Then either fi = s = 0, or 
[i = —s and k = 0. Either solution does not correspond 
to growing square-integrable modes. We can now assume 
u; > 0, and parameterise all solutions by s. We find 

fi(s) = vWw 2 , K{a) = f^p) ■ ( 21 °) 

We choose the principal branch of the square root, be- 
cause it maps Res > to Re [i > 0, that is, the growing 
modes are precisely the square-integrable modes. This 
choice also maps Res > to \k\ > 1. Therefore, we ex- 
clude all growing square-integrable modes if we restrict 
< 1 for all four Ki, or — 1 < Ki < 1 for real Ki. More 
details are given in Appendix iDl 

VI. CONCLUSIONS 

We have constructed families of generalisations of the 
BSSN and NOR variants of the ADM evolution equa- 
tions that are symmetric hy perb olic in the sense defined 
for second-order systems in |13j. This confirms the pre- 
vious result of |20| on the BSSN equations, without re- 
course to any first-order reduction, and generalises it by 
finding the most general energies for the main and con- 
straint system. This generalisation allows all characteris- 
tic modes to propagate with causal speeds, in particular 
with speeds (0, ±1) only. 

In our analysis of the BSSN equations we also clar- 
ify the role for hyperbolicity of imposing the algebraic 
constraints det 7^ = 1 and tr Ay = during the evo- 
lution. We find that the equations can be made sym- 
metric hyperbolic if these constraints are imposed con- 
tinuously, and strongly hyperbolic without imposing the 
constraints, but adding an extra term to the evolution 
equations. 

There is numerical evidence that densitising the lapse 
and imposing the trace constraint improves stability 



in moving single black hole simulations, even without 
imposing maximally dissipative constraint-preserving or 
boundary conditions 2j|. This is not surprising, as these 
changes make the evolution equations strongly hyper- 
bolic, and imposing the determinant constraint as well 
would make them symmetric hyperbolic. 

Our results go some way towards explaining why the 
BSSN system has been relatively successful in simulating 
black hole or neutron star binaries. It is possible that 
the Bona-Masso formulation j24|, a strongly hyperbolic 
first-order version of the Einstein equations that intro- 
duced variables similar to the T l or fi, has not been as 
successful because it is first order, which we expect makes 
it more susceptible to constraint-violating instabilities of 
the convergent type. 

The NOR system is basically the BSSN system without 
the conformal-traceless decomposition, and the similarity 
of our results for the two systems suggests that the NOR 
system shares all the advantages of the BSSN system, 
without the overhead of the extra variables K and 4> & n d 
extra constraints T = and D = 0. 

With symmetric hyperbolicity, we can make the initial- 
boundary value problem formally well-posed by imposing 
maximally dissipative boundary conditions. However, 
these boundary conditions are in general not compatible 
with the constraints, and so large constraint violations (of 
the convergent type) propagate in from the boundaries. 
This can be avoided by replacing some of the maximally 
dissipative boundary conditions on the main system by 
maximally dissipative boundary conditions on the con- 
straint system [l(| , and we have given details of how to 
do this for NOR and BSSN, for the case of a smooth 
boundary with tangential shift. Note that even when we 
have fixed the principal part of the field equations both 
the main and constraint energies still depend on a num- 
ber of free parameters, which appear explicitly in the 
boundary conditions. 

Except for two rather unphysical special cases, we have 
not proved that the initial-boundary value problem with 
constraint-preserving boundary conditions is well-posed. 
We have shown, however, that these boundary conditions 
rule out perturbation modes with unbounded growth, 
which is a key necessary condition for well-posedness £| . 
We plan to investigate a proof of well-posedness, and in 
parallel to examine stability in numerical experiments. 

Our discussion of both maximally dissipative and 
constraint-preserving boundary conditions assumes that 
the normal component (3 n of the shift vanishes at the 
boundary, as then the A = characteristic modes prop- 
agate along the boundary. This restriction allows for a 
shift that is everywhere tangential to a smooth boundary, 
and this could be used for example to employ corotating 
coordinates in the simulation of a binary system. The 
case of a general shift will be investigated in future work. 

All equations in this paper were derived using 
xTensor, an open-source Mathematica package for 
abstract tensor calculus, developed by JMM. It 
is available under the GNU Public Licence from 
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http : //metric . imaf f . csic . es/Martin-Garcia/xAct/ 
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APPENDIX A: THE KST FORMULATION 

The Kidder-Scheel-Teukolsky (KST) formulation is 
based on a reduction to first order of the ADM evolu- 
tion equations with a densitised lapse (DADM) with the 
auxiliary variables d k ij = dkjij- This gives rise to the 
auxiliary constraints Cijki = d^d^i — 0. The principal 
part of the evolution equations is 



dolij 



-2K, 



(Al) 



doKij ~ DADM + 77y -ff + C7 fc 'C fcfo)i , (A2) 
d d kij ~ DADM + r nk{i M j) +xiijM k . (A3) 

This system can be made strongly or symmetric hyper- 
bolic for certain ranges of the parameters a, 7, £, 77 and 
X- In particular, the Einstein-Christoffel (EC) system is 
the case 7 = 0, £ = —1,7? = 4,% = 0, densitizing the 
lapse with a = 1 in our notation. In this Appendix we 
want to point out that only a and 7 have counterparts 
in a second-order system. Q has a similar function to our 
parameter a, and 2r\ + x has a similar function to our 
parameter b, but these parameters vanish if we replace 
dkij by dhlij'- Cijki then vanishes identically, and dijk is 
no longer evolved explicitly by (| A3|> . but only implicitly 
by (HU). Comparing the KST system to NOR with the 
benefit of hindsight, one could say that the only indis- 
pensable effect of dkij is to introduce the divergence d k kj 
as an auxiliary variable. 



APPENDIX B: DIRICHLET BOUNDARY 
SYSTEM 

With 

— Kl = K2 = K3 = K4 = 1 

we have (in our example NOR system) 



(Bl) 



Fab = -2K AB , (B3) 
X = AK qq , (B4) 

X A = 2(f A - d n7A „) + (p - 2)d A (7nn + 7gg)-(B5) 

The boundary system 

d X = d A X A -2d A (d A j qq ), (B6) 

(B7) 

1 



d X A = ^d A X + 2d A F + 2d B F A B, 



d<){dAlqq) 



-d A X 



(B8) 



decouples from the bulk system. Note that all variables of 
the boundary system have parity +1 under the reflection 
rii — > — rii through the boundary. The boundary system 
is strongly hyperbolic with characteristic variables 

X ± = Ad m7qq -2X mT V6X, (B9) 
X a = X m + d mlqq , (BIO) 

which have speeds A = (± -y/3/2, 0) , as well as X p and 
dpj qq with zero speed. (It is also symmetric hyperbolic, 
but this is not required for well-posedness if the boundary 
is smooth without boundary, and therefore we do not give 
details here.) 



APPENDIX C: NEUMANN BOUNDARY 
SYSTEM 



With 

— Kl = K2 = K3 = K4 = — 1 (CI) 

we have 

F = 2/„ - d nlnn +{p- 2)d n { lnn + 7m ), (C2) 

U' Q = fn + ^d n ( 7nn + lqq ), (C3) 

Fab = d nlAB , (C4) 

-2d nlqq , (C5) 



X 

X A = ±K An . (C6) 
The autonomous boundary system is 

X = d A X A , (C7) 
X A - X -d A X + 2d B F AB + 20 A F 

-2d A % - 2d B d BlAn , (C8) 



9 (<9a7s™) = --d A X B , 
d U' Q = 0. 



(C9) 
(C10) 



Note that all variables of the boundary system have par- 
ity — 1. The boundary system is strongly hyperbolic with 
characteristic variables 

X± = 4d mlmn ~ X T V6X m , (Cll) 



F 



2K n 



(B2) 



Xq = dmlmn + T^X , 



(C12) 
(C13) 
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with speeds A = (±y/3/2, ±1, 0), as well as the zero speed 
variables Uq, d p -y mn and d p j pn . 

APPENDIX D: DETAILS OF THE MODE 
ANALYSIS 

An alternative approach to finding the range of real 
values of k that exclude growing modes would be to solve 
(I2()2|l and H2()9|) explicitly for the values of s and \i that 
are allowed for a given k, and then check if they corre- 
spond to a growing, square-integrable mode. For k = 
we have uj — and /i = — s. Therefore Re/i > implies 
Re s < 0, and vice versa, and so there are no growing 
square-integrable modes. For each n ^ we find two 
values of s and fi. For k > 0, we can parameterise them 
as 

D± : s = ±lo smhtp, ~oo < tp < oo, 

=>■ (i = ±u;coshtp, k = e~ 2ip . (Dl) 

The solution D~^_ with the upper sign and tp > has 
Re s > and Re /i > 0, and corresponds to growing square- 
integrable modes. To exclude these, we must demand 
K < 1. However, for k < there is a potential fallacy in 
this purely algebraic approach. 

Consider the complex s plane. For definiteness we 
place branch cuts from iio to ioo and from — iuj to — ioo. 
Consider the two contours C+ and C_ , parameterised by 
tp, which wrap around the upper and lower branch cut 
respectively: 

C± : s = ±iuj coship, — oo < tp < oo, 

=>■ fx = ±iu}s'mhip, K — —e 2ip . (D2) 



These give us the two possible values of fi and s for each 
real k < 0. From a point of view where we only con- 
sider real values of k, none of these modes are growing or 
square-integrable, and so we do not seem to obtain any 
restriction on negative n. However, we should consider 
modes with purely imaginary /i and s as limiting cases 
of growing square-integrable modes as Re s — > 0+ and 
Re /i — > 0+. Clearly this is the case for the modes lying 
on the contours C± (i.e. tp > 0) which are to the right of 
the branch cuts and therefore contiguous with Re s > 0, 
but not for the contours C± , which are to the left of the 
branch cut. For the modes on C±, \n\ > 1, and so they 
are suppressed by any boundary condition with |«| < 1. 
This is the range of K that we expect from the energy 
method. Another way of seeing that the modes on C± 
should not be considered is to note that our mode analy- 
sis is really derived from a Laplace transform in t so 
that we need a contour for the inverse Laplace transform 
that is to the right of all branch cuts (and singularities) 
in the s plane. 

The boundary of Re s > is given by the union of the 
three contours Cjj, Co and where 

„ . . 7T 7T 

C : s = lujsmip, -- < tp < -, 

/i = cj cos tp, k = e 2lip . (D3) 

This contour can also be used as a contour for the inverse 
Laplace transform. All contours in the s-plane discussed 
here and their images in the k plane are shown in Fig. ^ 
The shaded areas are Re s > and its image |«| > 1. Re s 
is bounded by by the union of the three contours C_ , Co 
and C.^ and Re \x by their images. This proves the claim 
made above that Res > is mapped to |k| > 1. 
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FIG. 1: Contours C and D in the complex s plane and their images in the k, plane under l|21()|l . Arrows point towards increasing 
parameter <p. The shaded area in the s plane is mapped to the shaded area in the k plane. 



